Chapter 6 Rasters in R: The terra package
For the longest time, when you wanted to manipulate raster files in R, the raster package was your tool of choice. And it still is a well-proven and tested alternative to the newer packages, one of which we will discuss today: terra.
6.1 Basic idea
terra implements two new object classes: SpatRaster and SpatVector. It is written and maintained by Robert Hijmans who also did so for the raster package. In that sense you can think of terra as the sequel to raster. Why did Robert Hijmans decide we need a new package? Well the problem of creating a widely used R package is that other packages start to use it and depend on it working they way it did when that package was written. So even if you don’t like certain aspects of you package anymore, you can’t just go ahead and rewrite them. Robert Hijmans thinks that raster has grown to be unnecessarily complex with over 200 functions. Its also slower than it could be and doesn’t support HDF5 files, which is a popular format for satellite data complex. So now terra is faster, simpler (functions have been streamlined) and more capable. However, when you know raster, you also know how to do most things in terra.
6.2 Creating a new Raster
The easiest way to specify a new raster is by calling rast() without any arguments.
## class : SpatRaster
## dimensions : 180, 360, 1 (nrow, ncol, nlyr)
## resolution : 1, 1 (x, y)
## extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (CRS84) (OGC:CRS84)
As you can see, the new object covers the whole earth and is projected using longitude and latitude on a WGS84 globe. The resolution is 1 by 1 which means each cell covers 1 degree of latitude and one degree of longitude. Of course, we can also specify the extend and resolution of the rasters we create. The following, raster for example, only covers the southern hemisphere.
## class : SpatRaster
## dimensions : 180, 360, 1 (nrow, ncol, nlyr)
## resolution : 1, 0.5 (x, y)
## extent : -180, 180, -90, 0 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (CRS84) (OGC:CRS84)
As you can see, this also automatically altered the resolution to 0.5. We cut the area in half so the resolution that corresponds to that direction is also cut in half. By explicitly setting the resolution, we can prevent this from happening.
Alternatively we can also specify the number of rows and cells we want our raster to have.
You can check the resolution of a raster with the res() function.
The function can also be used to compare the resolution of two (or more) rasters.
## [1] TRUE
We can also res() to change the resolution of a raster after creating or loading it.
## class : SpatRaster
## dimensions : 18, 36, 1 (nrow, ncol, nlyr)
## resolution : 10, 10 (x, y)
## extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (CRS84) (OGC:CRS84)
## class : SpatRaster
## dimensions : 2, 36, 1 (nrow, ncol, nlyr)
## resolution : 10, 100 (x, y)
## extent : -180, 180, -90, 110 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (CRS84) (OGC:CRS84)
The last thing we have not altered about this raster skeleton, is its coordinate reference system. If you only have an EPSG code for you desired CRS check out this website. It provides you with the corresponding PROJ.4 string.
## [1] "GEOGCRS[\"WGS 84 (CRS84)\",\n DATUM[\"World Geodetic System 1984\",\n ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n LENGTHUNIT[\"metre\",1]]],\n PRIMEM[\"Greenwich\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n CS[ellipsoidal,2],\n AXIS[\"geodetic longitude (Lon)\",east,\n ORDER[1],\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n AXIS[\"geodetic latitude (Lat)\",north,\n ORDER[2],\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n USAGE[\n SCOPE[\"unknown\"],\n AREA[\"World\"],\n BBOX[-90,-180,90,180]],\n ID[\"OGC\",\"CRS84\"]]"
# set to WGS 84
crs(x) <- "+proj=longlat +datum=WGS84 +no_defs "
# set to UTM 48
crs(x) <- "+proj=utm +zone=48 +datum=WGS84"
x## class : SpatRaster
## dimensions : 180, 360, 1 (nrow, ncol, nlyr)
## resolution : 1, 1 (x, y)
## extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=48 +datum=WGS84 +units=m +no_defs
6.3 Going beyond the hull
We can now create an empty raster anywhere on the world. However, as long as no values are in the cells this doesn’t get us anywhere. About time we introduce some values.
## [1] 100
## [1] FALSE
# now we fill the raster with increasing numbers, staring with 1
values(r) <- 1:ncell(r)
plot(r, main='Raster with 100 cells')
Let’s have a look at two actual raster file from my hard drive. You can download them here.
dem1 <- rast("data/DGM25_2905530.xyz")
dem2 <- rast("data/DGM25_2905540.xyz")
par(mfrow = c(1,2))
plot(dem1)
plot(dem2)
crs(dem1) = "+proj=utm +zone=32 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs "
crs(dem2) = "+proj=utm +zone=32 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs "As you can see they are two digital elevation models.
## class : SpatRaster
## dimensions : 161, 161, 1 (nrow, ncol, nlyr)
## resolution : 25, 25 (x, y)
## extent : 295987.5, 300012.5, 5535988, 5540013 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=32 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
## source : DGM25_2905530.xyz
## name : DGM25_2905530
## min value : 207.683
## max value : 511.909
6.4 Multiple layers in one SpatRaster
One SpatRaster can have multiple raster layers.
Think of the temperature in the same area in several years where each layer is one year.
Let’s simulate such an example.
# First we create a raster with temperatures between 10 and 15 degrees
r1 = rast(nrow = 25, ncol = 25)
values(r1) = runif(n = 25^2, min = 10, max = 15)
plot(r1)
We can combine the three rasters into one multi-layer object. With terra this is as easy as writing a simple vector
## class : SpatRaster
## dimensions : 25, 25, 3 (nrow, ncol, nlyr)
## resolution : 14.4, 7.2 (x, y)
## extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (CRS84) (OGC:CRS84)
## source(s) : memory
## names : lyr.1, lyr.1, lyr.1
## min values : 10.01700, 11.01700, 12.01700
## max values : 14.98155, 15.98155, 16.98155
Note that the class of a multi-layer object is the same as a single layer raster.
## [1] "SpatRaster"
## attr(,"package")
## [1] "terra"
Calling the generic plot function on such a multi-layer raster plots all layers, so be care full with this if you have many layers!

We can subset the object like a list.
## class : SpatRaster
## dimensions : 25, 25, 1 (nrow, ncol, nlyr)
## resolution : 14.4, 7.2 (x, y)
## extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (CRS84) (OGC:CRS84)
## source(s) : memory
## name : lyr.1
## min value : 10.01700
## max value : 14.98155
6.5 Raster algebra
SpatRaster objects can be passed to most algebraic operations in R such as + and - or sum() and abs().
You have just seen this above, when I added a degree of temperature to our simulated temperature rasters.
Below are some more examples
You can also add, substract, multiply and divide rasters with each other.
6.6 High-level functions
In most situations you will have rasters that are larger than what you need. You might be interested in the landcover of your study area but the data covers all of Europe. Conversely, you might also need to combine different rasters if the scale of your analysis is larger than the data you are using. Often this is the case when larger rasters are split in smaller areas to reduce their size.
In the first case you would use crop().
Provided with a raster and an extend, crop() will return only the part of the raster that lies within the extend.
An extend usually is a named vector with four elements: xmin, xmax, ymin and ymax.
It might for example be the bounding box you derived from the vector data you were analyzing in sf.
You can create this raster yourself if you know, or can derive, the relevant numbers.
An alternative is the drawExtent() function from the raster package.
This lets you click on an already plotted raster an returns the extend of the thus created rectangle.
Try the code below in your own R instance to see what I mean.


Nice!
Now have reduced the raster to the area we are actually interested in.
We can compute the mean hight or the sum of all hights with global().
## sum
## DGM25_2905530 100724.5
## mean
## DGM25_2905530 331.3307
We can also expand a raster with the extend() function.
That means we add rows and/ or columns with NA values.
The other way around we can also trim them, removing all cells with NAs in the margins.

Since all the NAs are just white, it’s hard to see what we actually did here.
Let’s give them some random value so we can see the cells we added.

We can now clearly see the added cells in brown.
With trim() we can remove them again.
na_cells = which(values(dem1.exp) == 300)
values(dem1.exp)[na_cells] = NA
dem1.trm = trim(dem1.exp)
plot(dem1.trm)
As mentioned before, sometimes you want to combine rasters.
For example our two DEMs are both part of a larger DEM that covers all of the German federal state Rhineland-Palatinate.
We can combine the two rasters with merge()

Whenever you want to use multiple raster in an operation, it is important that they have the same CRS and cell sizes.
You can increase cell sizes with aggregate() and reduce cell sizes with disaggregate().
In the first case, several cells are combined into one.
The number of cells is set with the fact argument.
When you supply one number it is used for both directions.
So if you call aggregate(x, fact = 2) and raster has the cell size 10m by 10m, the new raster will be 20m by 20m.
You can also provide fact with a vector that specifies the aggregation factor in each direction.
The third argument is fun which specifies how to calculate the new cell value.
If our new cell consists of four old cells (i.e., we used an aggregation factor of 2) we somehow need to combine those four numbers into one.
Common options here are the mean, the median, the modal, the maximum, the minimum, or the sum.
dem1_agg_2_mean = aggregate(dem1, fact=2, fun = "mean")
dem1_agg_3_mean = aggregate(dem1, fact=4, fun = "mean")
dem1_agg_4_mean = aggregate(dem1, fact=6, fun = "mean")
dem1_agg_5_mean = aggregate(dem1, fact=8, fun = "mean")
par(mfrow = c(2,2))
plot(dem1_agg_2_mean)
plot(dem1_agg_3_mean)
plot(dem1_agg_4_mean)
plot(dem1_agg_5_mean)
par(mfrow = c(2,2))
dem1_agg_mean = aggregate(dem1, fact = c(1,6), fun = "mean")
dem1_agg_min = aggregate(dem1, fact = 4, fun = "min")
dem1_agg_max = aggregate(dem1, fact = 4, fun = "max")
dem1_agg_sum = aggregate(dem1, fact = 4, fun = "sum")
plot(dem1_agg_mean, main = "directed aggregation")
plot(dem1_agg_min, main = "minimum")
plot(dem1_agg_max, main = "maximum")
plot(dem1_agg_sum , main = "sum")
When we want to go in the other direction, we disaggregate raster cells with the disagg() function.
Sadly, we cannot magically determine what values are the field truth for the new, smaller cells.
Instead each new cell gets the same value as it’s parent.

6.6.1 Extracting raster values with vectors
We often encounter situations where we want to extract the values of a raster that coincide with vector data we have. For example, we might have taken soil samples along a mountain and now we want to now the altitude of the samples or we are interested in the land cover of a river catchment.
To demonstrate this, we will need some vector data sets. For the point data, we can simulate points within the bounding box of the DEMs we used before.
# - extract extend of DEM
dem.ext <- ext(dem1)
# - convert extend to bounding box usable by sf
dem.bb <- st_bbox(dem.ext)
# - create 10 random points within bounding box
random_points <- st_sample(x = st_as_sfc(dem.bb), 10)
random_points <- st_as_sf(random_points, crs = 25832)
#- interactive map
tmap_mode("view")## tmap mode set to interactive viewing
Now we have ten points that are randomly distributed across the map.
If your repeat this with you own computer your ten points will be at other positions.
Now that we have, we can extract the raster values, i.e., the altitude, at the locations of our points.
This extraction is done with extract().
The first argument to extract() is the raster raster from which we want to extract values.
The second argument is the location of points, at which we want to extract the values.
## Warning: [extract] transforming vector data to the CRS of the raster
We get a non-spatial data.frame back that holds the row number of points and the elevation at these points. We can add the altitude values back to the points like this:
random_points$altitude <- ex.points$DGM25_2905530
tm_shape(random_points) + tm_dots(col = "altitude")Now we want to look at polygons instead of points. First, we need some polygons. Here, we use river catchments from Rhineland-Palatinate.
## Reading layer `catchments' from data source
## `C:\Users\jonat\Documents\001_Uni\002_teaching\online books\book_spatial_data_science_in_R\data\catchments.gpkg'
## using driver `GPKG'
## Simple feature collection with 2754 features and 55 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 4042743 ymin: 2874149 xmax: 4212789 ymax: 3094667
## Projected CRS: ETRS89-extended / LAEA Europe
Let’s have a look at the data.